import numpy as np
from scipy import interpolate #插值
import matplotlib.cm as cm #曲面的配色
import matplotlib.pyplot as plt #绘图
from mpl_toolkits.mplot3d import Axes3D #3D图

x0, y0 = np.mgrid[100:500:5j, 100:400:4j]#生成网格图5*4,5j代表5段，没有j代表间距为5
z0 = np.array([636, 697, 624, 478, 450,
               698, 712, 630, 478, 420,
               680, 674, 598, 412, 400,
               662, 626, 552, 334, 310]).reshape((4, 5))
z0=z0.T#注意z是5*4

f = interpolate.Rbf(x0, y0, z0)#由样本点生成三次样条插值
x1 = np.linspace(100, 500, 100)#np.linspace(start, stop, num=50)
y1 = np.linspace(100, 400, 100)
z1 = f(x1, y1)#插值结果

x1, y1 = np.meshgrid(x1, y1)
fig = plt.figure()
ax=Axes3D(fig)#3D对象
#ax=fig.add_subplot(111,projection='3d')
#ax= plt.subplot(1, 2, 1, projection='3d')
surf = ax.plot_surface(x1, y1, z1, rstride=1, cstride=1,
                         cmap=cm.coolwarm, linewidth=0.5, antialiased=True)#rstride行之间的跨度，cstride列之间的跨度cmap颜色映射表，默认是“rainbow”
plt.title('Interpolation-2D(The peak: {:3f})'.format(np.max(z1)))
ax.scatter(x0,y0,z0,c='r')#三维散点图
ax.contour(x1, y1, z1,zdir='z',offset=300)#投影到z平面的z=300上
ax.set_zlabel('Z')  # 坐标轴
ax.set_ylabel('Y')
ax.set_xlabel('X')
plt.colorbar(surf, shrink=0.5, aspect=5)
plt.show()

#最高点处的地址
x1[np.where(z1==np.max(z1))] #array([144.44444444])
y1[np.where(z1==np.max(z1))] #array([200.])